function varargout = flowc_from_imgs(VBas, imgs, dx, dy, dt)
%FLOWC_FROM_IMGS Estimation of flow coefficients from image sequence
%
%   alpha = FLOWC_FROM_IMGS(VBas, imgs, dx, dy, dt);
%
%       Estimates the flow coefficients w.r.t to a given flow basis,
%       from an image sequence. 
%
%       The estimation is formulated into an optimization problem as
%       follows:
%
%           minimize (1/2) * || predicted changes - observed changes ||^2.
%
%       Input arguments:
%       - VBas:         The velocity bases of the flow family
%       - imgs:         The sequence of images (size: [h, w, n])
%       - dx:           The interval between neighbor pixels (along x)
%       - dy:           The interval between neighbor pixels (along y)
%       - dt:           The time interval between consecutive frames
%          
%       Output argument:
%       - alpha:        The estimated coefficient vector. [K x 1].
%
%   [H, f] = FLOWC_FROM_IMGS( ... );
%
%       Returns the accumulated statistics. 
%
%       - H:    The accumulated precision matrix
%       - f:    The accumulated potential vector
%
%       They can be used in Bayesian estimation. For MLE, the optimal
%       coefficient is given by H \ f.
%

% Created by Dahua Lin, on April 5, 2012
%

%% verify input arguments

if ~is_vbasis(VBas)
    error('flowc_from_imgs:invalidarg', ...
        'VBas should be a velocity field basis.');
end

h = VBas.size(1);
w = VBas.size(2);
K = VBas.K;

if ~(isfloat(imgs) && isreal(imgs) && ndims(imgs) == 3 && ...
        size(imgs,1) == h && size(imgs,2) == w)
    error('flowc_from_imgs:invalidarg', ...
        'imgs should be a 3D array with size(imgs,1) == h and size(imgs,2) == w.');
end
n = size(imgs, 3);

if ~(isfloat(dt) && isreal(dt) && isscalar(dt))
    error('flowc_from_imgs:invalidarg', 'dt should be a real scalar.');
end

%% main

H = zeros(K, K);
f = zeros(K, 1);

for i = 1 : n-1
    
    I_cur = imgs(:,:,i);
    I_next = imgs(:,:,i+1);
    
    [Gx, Gy] = imgrad(I_cur);
    Gx = Gx * (1/dx);
    Gy = Gy * (1/dy);
    
    I_cur = I_cur(:);
    I_next = I_next(:);
    
    B = calc_base_changes(VBas, Gx, Gy);
    c = (I_next - I_cur) *(1/dt);
    
    H = H + B' * B;
    f = f + B' * c;    
end

if nargout <= 1
    alpha = H \ f;
    varargout = {alpha};
else
    varargout = {H, f};
end


%% auxiliary functions

function B = calc_base_changes(VBas, Gx, Gy)

Cx = bsxfun(@times, VBas.Bx, Gx(:));
Cy = bsxfun(@times, VBas.By, Gy(:));
B = - (Cx + Cy);







